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We investigate the static and dynamical behavior of ID interacting fermions in disordered Hub- 
bard chains, contacted to semi-infinite leads. The chains are described via the repulsive Anderson- 
Hubbard Hamiltonian, using static and time-dependent lattice density-functional theory. The dy- 
namical behavior of our quantum transport system is performed via an integration scheme available 
in the literature, which we modify via the recursive Lanczos method, to increase its efficiency. To 
quantify the degree of localization duo to disorder and interactions, wo adapt the definition of the 
inverse participation ratio to obtain an indicator which is both suitable for quantum transport ge- 
ometries and which can bo obtained within density-functional theory. Lattice density functional 
theories are reviewed and, for contacted chains, we analyze the merits and limits of the coherent- 
potential approximation in describing the spectral properties, with interactions included via lattice 
density functional theory. Our approach appears to able to capture complex features due to the 
competition between disorder and interactions. Specifically, we find a dynamical enhancement of 
delocalizatiou in presence of a finite bias, and an increase of the steady-state current induced by 
inter-particle interactions. This behavior is corroborated by results for the time-dependent densities 
and for the inverse participation ratio. Using short isolated chains with interaction and disorder, 
a brief comparative analysis between time-dependent density-functional theory and exact results is 
then given, followed by general conclusive remarks. 

PACS numbers: 31.15.ee, 72.15.Rn, 72.10.Bg, 71.10.Fd 



I. INTRODUCTION 

In many physical phenomena, practical limitations hin- 
der a complete knowledge of all the degrees of freedom 
involved. Nanoscience has adopted such apparent short- 
coming as its central paradigm, by exploiting the notion 
of a small system coupled to a macroscopic environment. 
A case in point is represented by nanoscale transport phe- 
nomena, where two (or more) macroscopic leads are con- 
nected to small central devices (quantum constriction)^'^. 

Such devices, whose size ranges from that of few atoms 
(as in short nanowires or small molecules) to that of sev- 
eral repeated large molecular units, attract scientific in- 
terest because they are seen as possible candidates for 
novel electronic, spintronic, or qiiantum computation de- 
vices, to mention a few"^. This, in turn, requires a thor- 
ough understanding and control of the decoherence pro- 
cesses which can affect carrier propagation and manipu- 
lation in the device region. 

In this work we consider two of such processes, namely 
disorder and inter-particle interactions (thus leaving out 
other important decoherence mechanisms, e.g., lattice vi- 
brations). How interactions and disorder affect the con- 
duction properties of materials has been intensively inves- 
tigated over several decades'*"^*', and significant progress 
has been made. However, some issues remain at a consid- 
erable extent open, e.g. the real-time dynamics of sam- 
ples with disorder and interactions. 

Starting with the seminal paper by Anderson^, lattice 
models have had an eminent place in the study of dis- 
ordered systems with and without interactions. While 
a large fraction of the literature on disordered inter- 



acting lattice models concentrates on the equilibrium 
regime (for both finite and extended systems), more re- 
cently the time-dependent properties have also been ex- 
amined, primarily for finite samples^^'^^. In between 
the finite/infinite-system categories, a third one is rep- 
resented by small disordered samples connected to semi- 
infinite homogenous reservoirs^®, of relevance to quantum 
transport phenomena. 

This paper looks into some aspects of the transport 
properties of ID interacting fermions in disordered lattice 
systems, using static^" and time-dependent density func- 
tional theory^i (DFT and TDDFT, respectively). Static 
and time-dependent DFT are in principle exact reformu- 
lations of the (time-dependent) many-body problem^^, 
where the key variable is the one-particle density n, and 
a central ingredient is the exchange-correlation (XC) po- 
tential Vxc (recent comprehensive reviews of the subject, 
are Refs. 23-25). The XC potential embodies the com- 
plexities of the many-body problem. In this contracted 
description, 

Vxc ^ 

highly non-trivial functional of the 
density (in TDDFT, where time enters explicitly into the 
formulation, such functional dependence includes the en- 
tire history of the density n, i.e. memory effects). In 
general, the exact Vxc is not known, and approxima- 
tions are introduced. A simple but not always adequate 
prescription is the so-called (adiabatic) local density ap- 
proximation, where the XC potential depends only on 
the local (time-dependent) density. This amounts to ne- 
glecting non-locality in space (and memory effects in the 
TD case) in Vxc- As a result of this oversimplification, 
in some practical applications an accurate description of 
dynamical inter-particle correlations may be lacking. 
The application of static density-functional theories 
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to lattice models started almost thirty years ago^^ 
and in the last decade this approach has been further 
developed^""^^. The use of lattice TDDFT to describe 
the non-equilibrium dynamics of Hubbard-like models is 
a rather new topic^'*"^'', and some of its formal aspects 
are still under scrutiny (firm conceptual ground has been 
established for related lattice approaches which use the 
lattice bond-current as the basic variable^^"*^). However, 
there is a growing body of evidence (see Sect. Ill) sug- 
gesting that the lack of a rigorous formulation (an issue 
which is likely to be resolved in future) could be of no 
significant practical consequence. 

After the above, somewhat lengthy, considerations, we 
can now define the motivations behind this work and the 
chosen methodology. Our focus is on finite chains con- 
tacted to semi-infinite leads, with short-range interac- 
tions and disorder present only in the chains (the "de- 
vice"). Even within these narrow boundaries, the issues 
which can potentially be addressed are many, but we will 
only touch upon a few of them, and going in no great de- 
tail. In this respect, our work is somewhat exploratory 
in character, since we also describe some methodologi- 
cal developments, that we found necessary when using 
(TD)DFT for disordered lattice systems in the quantum 
transport regime. Concerning the chosen methodology, 
we note that, compared to other lattice approaches, lat- 
tice (TD)DFT is well suited for any dimensionality and 
is relatively inexpensive from the computational point of 
view (since it deals with single-particle orbitals). These 
are attractive features when, for example, one needs to 
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Eq.(l) describes a central chain of length L (the lattice 
sites with \ <l <V) connected to a left and a right ID 
lead [sites with I < 1 and I > L, respectively. The third 
and fourth term in Eq. (1) represent the time-dependent 
bias in the leads (r is the time variable), which is applied 
at time r > [often, in the literature, the leads are also 
referred to as the source (S) and drain (D), hence the 
subscripts S, D in the bias terms in Eq.(l)]. For the con- 
tacted chain, the hopping term Vij+i = V always, i.e. 
we employ transparent boundary conditions (hereafter, 
V = 1 is taken as the energy unit). The Hamiltonian 
of the isolated chain is obtained from the general one 
by putting Vb,i = Vl,l+i = in in Eq.(l), and retain- 
ing only the sites labeled hy 1 < I < L. Looking more 
closely to the chain part of the Hamiltonian, we have 
Hubbard-like interactions (the term proportional to U ; 
we set Z7 > 0) and time-dependent onsite energies wi{t), 
which is convenient to separate into a static and time- 
dependent part: wi{t) = e; +vi{t). In the presence of 



perform configuration averages of time-evolved quanti- 
ties in the long-time limit (such as when approaching the 
steady state regime in a quantum transport setup with 
disordered samples). 

The outline of the paper is as follows: In Sect. II, we 
describe the lattice model system that we employ in our 
work. In Sect. Ill, we present the lattice (TD)DFT 
formalism, and describe how to obtain the XC poten- 
tial from the exact solution of the ID Hubbard model. 
This is accompanied by a review of the inherent litera- 
ture, to illustrate the developments and applications oc- 
curred so far within this approach. In Sect. IV, which 
deals with disorder, we first discuss the inverse participa- 
tion ratio, then we introduce a formulation for contacted 
chains based on the Coherent-Potential Approximation 
and DFT. In Sect.V we start by briefly reviewing lat- 
tice TDDFT approaches to quantum transport. Then, 
we present in some detail a method recently proposed in 
the literature^^ , followed by a description of our modifi- 
cations to it, to increase its efficiency. Some technical de- 
tails relative to Sects. IV,V are relegated to Appendixes 
I and II. In Sect. VI, we report and discuss our results, 
for static and non-equilibrium regimes. Our conclusions 
are in Sect. VII. 

II. THE MODEL 



In standard notation, the lattice systems considered in 
this paper are described by the following Hamiltonian: 



nia + bsir) ^ nia + bniT) ^ nia- (1) 

i<l;<T l>L;a 
I 

disorder, the e;:s are distributed according to some disor- 
der probability distribution. In this work, we use primar- 
ily the box disorder distribution, i.e. e; e [—W/2,W/2], 
but, we will sometimes consider binary disorder, where 
e; = ±W/2. In both cases, W fixes the strength of the 
disorder. The chain Hamiltonian is a finite-size realiza- 
tion of the so-called Anderson-Hubbard model (AHM)^°, 
one of the most used models to study strongly corre- 
lated and disordered systems^^. The AHM generalizes 
the standard Hubbard Hamiltonian^^ to inhomogeneous 
(and, in our case, possibly time-dependent) situations. 
That is, V(t) = X^icr '''(''') "'t describes a local (in space 
and time), time-dependent potential in the chain. In the 
static case (i.e. before the systems starts to time-evolve), 
all vi{t) = 0. Furthermore, the usual Hubbard model for 
the chain**^ is recovered when Wi{t) = 0, whilst, when 
f/ = but e; ^ 0, the chain is described the so-called 
Anderson model of disorder^. 
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III. (TD)DFT FOR LATTICE MODELS 
A. General aspects of lattice (TD)DFT 



can write for the ground-state total energy^*'^^: 
E[n, Vext] = To[n] + Enin] + E^cin] + ^ Vext{i)nu (2) 



A DFT based on the site occupation numbers nn was 
introduced more than two decades ago, to describe some 
ground state properties of the Hubbard modcP^"^*. An 
exact LDA (based on the Bethe-Ansatz) for the inhomo- 
geneous ID Hubbard model was first considered in Ref. 
29. Further significant progress came when an explicit 
and simple expression for the XC functional based on the 
Bethe-Ansatz was provided^*^, and practically used to in- 
vestigate different inhomogeneous Hubbard-type models. 
In subsequent work, The LDA based on the Bethe-Ansatz 
for Vxc was scrutinized against exact results^^^^^''*^''''*, 
providing energies, particle densities and entropies with 
an accuracy within a few percents. 

Recently, lattice DFT has also been used to determine 
the polarizability of the ID Hubbard model^^, and also to 
study the entanglement entropy of the Hubbard model^*. 
Furthermore, explicit analytical expressions for the XC 
potentials in small clusters can be found in Ref. 46, while 
the role of the temperature on v^^ has been discussed in 
Ref. 47. Lattice (TD)DFT has also been used to inves- 
tigate ultracold atoms on an ID optical lattice^^'^'"'^^"^^. 
These systems permit to study different ground-state and 
non-equilibrium scenarios for the Hubbard model^^ with 
high accuracy (because a precise tunability of the lattice 
parameters is possible), more directly and easily than in 
solid-state experiments. In 2D, the Hubbard model has 
been investigated via DFT on the graphene lattice^^. To 
date, the simple cubic lattice is the only 3D case con- 
sidered in the literaturc^^, with the ground-state energy 
of the uniform system computed within dynamical mean 
field theory (DMFT)^'*'^^ Differently from the ID case, 
here a discontinuity in Vxc appears only for U > U^°**, a 
DFT description of the onset of the Mott-Hubbard metal- 
insulator regime at a finite U. 

While lattice DFT rests on rigorous grounds^^, there 
is at present no direct formulation of the Rungc-Gross 
theorem for lattice TDDFT, as discussed in Ref. 38; 
early discussions of t;o-representability on the lattice can 
be found in Rcfs. 34-36. On the other hand, with the 
bond-current as the basic variable, a rigorous formulation 
on the lattice becomes possible^^"^*^. Finally, it should 
also be noted that in ID systems, as those considered in 
this paper, there is a one-to-one correspondence between 
densities and currents, and thus TDDFT rests on solid 
grounds. Lattice (TD) DFT has also been considered in 
context of work on quantum transport geometries; this 
aspect will be examined further below. 



B. Formulation 

In this paper, we confine ourselves to the non magnetic 

ID case; we review here the actual formulation for spin- 
independent (TD)DFT. In standard DFT notation, we 



where v^xtii) = is the static external field, and 
TqIu] and Eh = j'^iUiu'l arc, respectively, the non- 
interacting kinetic energy and the Hartree energy, with 
rii ~ rii^. To perform a local density approximation, 
Exc is obtained from a homogeneous reference system 
(Hubbard model): 



Exr. = E — Tc\ — E 



H 



(3) 



To obtain Vxc, one takes the derivative of the XC energy 

per site Cxc = Exc/L with respect to the density (in the 
general case, this should be a functional derivative): 



de xc{n, U) 
dn 



(4) 



For bipartite lattices, Cxdn, U) = Cxd^ — n, U) in the 
entire density range [0, 2] and thus Vxc{n) = —VxcC^ — n). 
Finally, a local-density approximation is defined: 



S) = Vxc{ni). 



(5) 



In ground-state DFT-LDA calculations, the XC potential 
obtained in this way is used to to solve self-consistently 
the Kohn-Sham (KS) equations 



{i + VKs)'PK = Sk^k 



(6) 



where i denotes the matrix for the single-particle hop- 
pings between nearest-neighboring sites, and ^p^ is the k- 

th single-particle KS orbital, with rii = X^kgocc Iv'k(OP- 
The effective potential matrix is diagonal: {vKs)ii = 

VKs{i) = VH{i) + Vxcii) + Vext{i), with VH{i) = lUiTli 

being the Hartree potential. 

If DFT can be a viable route to describe the ground- 
state properties of Hubbard-type models, then the lattice 
Kohn-Sham (KS) equations could be propagated in time, 
to get a TDDFT description of the dynamics of lattice 
systems. For ID Hubbard-type Hamiltonians, work in 
this direction was performed in Refs. 56-58 for the linear 
response regime. A TDDFT approach to the real-time 
dynamics of the Hubbard model out of equilibrium was 
first considered in Ref. 35, where exact results for the 
density and the XC potentials were compared to those 
obtained by solving the time-dependent KS equations 



{t + vks{t)) <Pk{t) = idrif^ir) 



(7) 



In general, VKsihr) = VH{i,T) +Vxc{i,T) +VextihT) de- 
pends non-locally on the density via Vxc- The adiabatic 
local density approximation (ALDA)^^ to the XC poten- 
tial is then obtained with the prescription v^^^^{i, r) = 



.,LDA 



(ni(r)), where the TD density is given by ni(r) = 
J2"tt^ !</'«(«, t)!^- An ALDA for the Hubbard model was 
first introduced in Ref. 35, with the treatment limited 
to spin-compensated systems, while the spin-dependent 
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case was presented in Ref. 37, where TDDFT results 
and time-dependent DMRG (tDMRG) results were com- 
pared. 

For finite systems, a study more focussed on the role of 
non-local and memory effects beyond the ALDA was per- 
formed in Ref. 60, via the ALDA, exact, and Kadanoff- 
Baym time-evolution in small cubic Hubbard clusters. 
The Kadanoff-Baym equations (KBE), with a many- 
body perturbation-theory approach to the self-energy, 
permit to take into account non-locality and memory ef- 
fects on equal footing. Such comparisons showed that an 
ALDA coming from the appropriate (strongly correlated) 
reference system can perform well in many instances, (es- 
pecially for slow perturbations) but, quite generally, it 
will fail for fast perturbations, or very strong interac- 
tions. 

We conclude this Section with a remark about no- 
tation: the one adopted throughout the paper is fully 
consistent with the continuum case, i.e. as if the Hub- 
bard interaction was treated as spin-independent: U = 
\YliiUi{n^ — fii). However, when the interaction is 
rewritten as in Eq. (1) the interaction is effectively kept 
among opposite spins (i.e., treated as spin-dependent)^^ 



and the exchange has been removed at the Hamiltonian 
level. Thus, e^c and v^c in Eqs. (3,4) contain only corre- 
lation, and the Hartrec plus exchange potential is Uni/2. 
More aptly, our DFT quantities could have been called Co 
and Vc but, following a common practice in the literature 
on lattice (TD)DFT, we still denote them by e^c and Vxc- 



C. Obtaining Vxc for the ID Hubbard model 

According to Eq.(3), to construct a LDA in ID 
we need the exact ground-state energy of the infinite 
homogeneous ID Hubbard model, where the hopping 
= V and the interaction is present at all sites. 
This requires^^'^*^'^^'"*^ to solve the coupled Bethe-Ansatz 
equations for the charge and spin distribution functions 
[p{x) and cr{x), respectively]^^. An ALDA is then easily 
obtained^^'^'^, making Vxc to become a function of the 
instantaneous local density. In the non-magnetic case 
considered here (where the spin-up and spin-down den- 
sities are equal, i.e. n-^ = = n/2), the Bethe-Ansatz 
equations read 



, , 1 cosx /■"'"°° u/4 , , , 

"^^^ = ^ L (V4)^ + (.-sin.)^^(^)^^ - ^ L iu/2r + iy-y'f^' ' 



with u = U/V. The functions p{x) and (t(x) are related 
to the charge n = n^ + ni and spin-down densities via 

/+Q r+oo 
p{x)dx = n, / (T{y)dy = n^, (10) 
-Q J-oo 

from which the integration limit Q is determined. The 
ground state energy density for n^. = is given by 

r+Q 

e{n^ = ni,U) = —2V / p{x)cosxdx. (11) 
J-Q 

Eqs. (8-11) are the prescription used in this work to deter- 
mine Cxc- The numerical solution of Eqs. (8,9) was ob- 
tained via a self-consistent procedure, with Q adjusted at 
each iteration via the normalization condition in Eq.(lG). 
Numerical integrations were performed with a 128-point 
Gauss-Legendre quadrature and, for each U, exc{n, U) 
was obtained at the nodes of a uniform mesh for the den- 
sity n. To obtain Vxc at each node of the density mesh, we 
computed Sexd^n with a 5-point numerical derivative. 
To calculate Vxc at off-node densities, a linear interpola- 
tion between the the closest nodes was employed. The Vxc 
thus obtained is discontinuous at half-filling, as it should 



be for the homogeneous ID Hubbard model; however, for 
a finite interacting system contacted to non-interacting 
leads, the discontinuity of the exact v^c becomes slightly 
smoothened (this was indicated in Ref. 63, using support 
from small Anderson clusters, and fully discussed in Refs. 
47. G4 66). According to these considerations, and also 
for numerical convenience, the XC potential was slightly 
smoothened in our actual calculations. 

IV. THE ROLE OF DISORDER 

Methodologically, the way we numerically deal with 
disorder effects in this paper is straightforward, since in 
most cases we limit our analysis to the arithmetic (con- 
figuration) average of specific quantities: the inverse par- 
ticipation ratio (IPR), the density and the current den- 
sity. In most cases, the numerical configuration averages 
are performed over an incomplete set of configurations. 
Sects. IV A and IV B constitute an exception. In these 
sections, which deal with ground-state quantities, and for 
the case of binary disorder, we also perform complete nu- 
merical averaging over all the configurations, to provide 
benchmark results. For brevity, in the following, com- 
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plctc numerical averaging will be referred to as 'exact 
averaging'. 

A. The inverse participation ratio 

A quantity often used as an indicator of localization 
in a system is the so-called inverse participation ratio 
(. The original definition^'' of introduced for non- 
interacting disordered systems, characterizes a given one- 
particle quantum state as follows: 



M 



M 



(12) 



where = is the density at site i and the sums 
extend to all the sites M in the system. For completely 
localized states (when 7^ at only one site) wc get 
Co = 1, while Co is smallest for delocalized states. To 
deal with interacting systems, suitable modifications of 
Eq.(12) are, for example. 



M 



M 



EAn?/(EAn.)^ 



M 



M 



C2H = E"?H/[E"^H]'- 



(13) 



(14) 



The use of Ci is convenient when dealing with small sys- 
tems with discrete many-body levels^^. In this case, for 
N particles, Aui is the difference between the ground- 
state densities with + 1 and N particles, a clear opera- 
tional prescription for finite systems. Using (^2^^ amounts 
to consider the density of states as obtained from the 
one-particle propagator, since ni(w) = —QGii (a;)/7r. It 
should be noted that most investigations of the IPR are 
done numerically, for finite systems. Using C2 requires in- 
troducing a finite artificial broadening 7, and employing 
a finite-size scaling analysis, to assess the role of -y'^o^''^. 

In this work, we study finite disordered systems (short 
chains) contacted to semi-infinite homogenous leads. 
While the definition of IPR via Eq.(14) is in principle 
suitable, for our lead-device-lead system we are faced 
with two issues: i) a sum is implied over all the sites 
M in the system, including the leads; this considerably 
increases the difficulty of the calculation ii) more fimda- 
mentally, the IPR defined in this way can reflect the prop- 
erties of the leads rather than the actual systems, since 
the lead-contribution can dominate the sums in Eq.(14). 
In view of this, we suggest the following possible alterna- 
tive definition of the energy dependent IPR: 



C3 = E"?H/[E"*('^) 

iec iec 



(15) 



i.e. the sum is confined to the region of the device. This 

modified definition of the IPR will be further analyzed 
and compared to the standard one in Sect. VIA. 



B. The Coherent Potential Approximation 

Among the possible theoretical approaches to disorder, 
an important place is occupied by the Coherent- Potential 
Approximation'''*"^^ (CPA), which introduces a simpli- 
fied, approximate treatment of disorder averaging'^^"''^. 
A major appeal of CPA in its simplest formulation is 
the pedagogical value, a relative analytical simplicity and 
ease of numerical implementation, together with the abil- 
ity to give results for several quantities (e.g. ground state 
energies, transport properties, phase stability, photoe- 
mission) which are generally in broad agreement with 
experiment. Furthermore, the CPA becomes exact in 
the D = 00 limit ''^; in finite dimensions, it has been nu- 
merically tested against accurate mimcrical averaging^" 
and it has also been used in conjunction with many- 
body perturbation theory (see e.g. Ref. 81). It has 
also been shown*^ that, when possible, numerical aver- 
aging based on small optimized supercells can give re- 
sults considerably superior to the CPA. Furthermore, it 
should be pointed out that the CPA shows significant 
limitations when describing quantities where spatial cor- 
relations among different sites cannot be neglected. A 
discussion of the properties of CPA (its limits of valid- 
ity, extensions, applications, etc.) is outside the scope of 
this work, and here we provide a short derivation which 
combines CPA and lattice DFT for quantum transport 
geometries. To this end, we must slightly modify the 
standard treatment '^^"''^ , and adapt it to the case of a 
finite disordered sample in the presence of homogeneous 
semi- infinite contacts, and within a lattice DFT context. 
We consider the case of diagonal disorder, and special- 
ize to a finite sample (chain of length L) of a random 
binary alloy, with species A and B and concentrations 
ca = Na/L and cb = ^ — ca, respectively. In the 
quantum transport geometry we study, the chains are 
connected to non-interacting leads, as in Scct.V. For 
such chains, complete disorder averaging for a given con- 
centration requires [j^^) configurations, and already for 
L = 14, 15 (as considered in this work), this number is 
rather large. In fact, performing time-dependent quan- 
tum transport calculations based on complete numerical 
averaging is computationally very demanding, and one 
is bound to use a much reduced (and incomplete) nu- 
merical sampling. This latter strategy is the one mainly 
adopted in the paper. To assess the scope of CPA,we 
limit ourselves to the ground state, when no bias is ap- 
plied. Numerical results relative to this formulation are 
presented in Section VI B. 

In matrix notation (in the site indexes) , and in absence 
of disorder and magnetic effects, the retarded one-body 
Kohn-Sham propagator for a chain connected to a left 
and a right lead can be written as 



(16) 



wl - Hks - ctl(w) - (7^(0;) 

where the matrix indexes of g^s label the sites of the 
chain, is the self-energy operator from the left 
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(right) lead*^, and Hks = t + vks, accounts for the 
kinetic and potential Kohn-Sham (KS) operators. In the 
presence of binary diagonal disorder H Hks + V, 
where V = X^i^j^i- Here, hi = Y.a^^'y- ^or e^, the 
probability distribution is P{ti) = CA6{ei — ca) + (1 — 
c-Aj^i^i — es)- The CPA provides a prescription to deter- 
mine the effect oiV. In an exact fashion, we can equiv- 
alently write for the propagator (Gks) averaged over all 
disorder configurations (the local dependence on co is not 
shown): 

{Gks) = 9KS + 9KS S {Gks) (17) 
{Gks) = gKS + gxs (T) 9ks, (18) 

(T is the T-matrix of the potential.) Inserting Eq.(18) 
in Eq.(17), and after some simple manipulations, we get 



This equation must be solved for each site in the chain, 
once the local propagator {GKs)ii is known. The latter 
is in turn determined from Eq.(16), after the replacement 

Hks ~^ Hks + S"-^^"^ is made, and after the dependence 
of Hks on the densities has been taken into account via 



rii 



r 

J —CSC 



9(GKs(w + iO+)), 



(24) 



with /i the chemical potential (here, as in the rest of 
the paper, we work at zero temperature). While it is 
certainly possible (and often necessary) to improve over 
the CPA®^, in this work we aim at qualitative insight, 
and in Sect VI B we present results obtained with the 
simple local formulation of CPA and the self-consistent 
set of equations Eqs.(16,22,24)^*5. 



S = (r) l+5ifs(T) 



(19) 



In the CPA, the correlations among different scatterers 
ar{^ taken into accoimt by assuming an effective medium 
for which the single site scattering is zero on average. To 
see how, we start with a specific disorder configuration V , 
for which Gks = 9ks + 9ksVGks- Then, we subtract 
from both sides the quantity f^GKS (with S yet to be 
specified). This gives {g^s — S) = 1 + (V^ — Ti)Gks and, 
if we choose (^^s ~'^) — {Gks)~^, i-e. as in Eq.(17), we 
finally get 



Gks = {Gks) + {Gks){V - T,)Gks 
= {Gks) + {Gks) (t) {Gks) 



(20) 



Performing the average of Eq.(20) over different configu- 
rations , we note that it must be (T) = 0, the key exact 
condition for the T-matrix. 

The CPA makes two assumptions: i) E is diagonal in 
the site- indexes, ^fj^^ = ^ij^f^^^ and so is the per- 
turbation V — T,; ii) instead of (T) = 0, one imposes 
a simpler approximate constraint, i.e. that the average 
local T-matrix at the i-th site is zero: 



{t{i)} = = CA{tA{i)) + CB(tB(«)), 



where 



tA(B)(«) 



£A{B) 



^CPA 



l-(Gx5)i4eA(B)-Sr^] 



(21) 



(22) 



In contrast to the usual treatments here all quanti- 
ties in Eq.(22) depend on the site- index, since our system 
exhibits no disorder in the leads. Very recently, and in- 
dependently, a similar formulation has been provided in 
Ref. 84. 

Inserting Eq.(22) in Eq.(21), and performing simple 
algebra, we arrive to an equation for S?-^^: 



^CPA ^ 
- ^A 



(1 - CA){eB - ca) 
l-(Gxs)u[eB-Sr^] 



(23) 



V. TIME DEPENDENT QUANTUM 
TRANSPORT (TDQT) AND TDDFT 

Theoretical approaches to quantum transport can be 
broadly grouped according to different criteria, e.g. if 
they are based on a steady-state or time-dependent 
formulations, if they use ah initio or model Hamilto- 
nian methods, or according to which mathematical tech- 
nique is employed: non-equilibrium-propagator, linear- 
response, wavefunction-scattering, etc. Here we consider 
time dependent quantum transport (TDQT), which per- 
mits to follow the system during its time evolution after 
a bias has been applied. In this way, steady-state, tran- 
sient and a.c. currents can all be considered on equal 
footing and, in the presence of dissipation, history de- 
pendence (memory effects) are also accounted for. A vi- 
able strategy to TDQT is to consider large but finite 
systems^''. Via an initial spatial imbalance of particles, a 
quasi-steady state current can be established. Recently 
this approach has also been used to describe bosonic 
and fermionic transport of ultra-cold atoms in ID optical 
lattices^S'S^. 

A different formulation, the one used here, consid- 
ers a central region initially connected to semi-infinite 
leads^°~^^. This "contacted" approach has been used to 
introduce a TDDFT description of TDQT''l''^'^ and the 
practical applicability of a TDDFT scheme has also been 
shown^^. Furthermore, classical nuclear degrees of free- 
dom have also been included in the approach^''. 

In a TDDFT approach to TDQT^^'^^^ a key quantity 
is the XC potential. In Ref. 95, a combination of DMRG 
and lattice DFT was used to gain insight into the ex- 
act ground-state XC functionals for a correlated-electron 
model system coupled to external reservoirs. A compari- 
son of lattice DFT and DMRG in transport has also been 
provided in Ref. 96, whilst a study of the role of spin in 
the XC potential can be found in Ref. 97 (for treatments 
based on ground-state current DFT for lattice models, 
see Refs. 98,99). 

The effect of a discontinuity in Vxc in a TDDFT de- 
scription of TDQT was examined within lattice TDDFT 
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in the ALDA approximation^^. Following the time- 
evolution of a single Anderson impurity attached to two 
biased leads, a dynamical notion of the Coulomb blockade 
was then presented. This emerges also from a description 
based on time-dependent, bond-current DFT^^. In Ref. 
63, it was also pointed out that for a single Anderson 
impurity, the exact v^c is a sharp (at half-filling ) but 
smooth function of the density. 

Subsequently, a comparison between ALDA, tDMRG 
and KBE for TDQT in lattice systems was presented in 
Ref. 100, showing that the ALDA can give accurate den- 
sities but overestimated currents, due to the neglect of 
non-local effects in the leads. We finally mention that, 
very recently, different research groups^^'^^'^*^^ indepen- 
dently pointed out that suitable XC potentials permit 
a (TD)DFT description of the Kondo effect, and also 
examined in detail the broadening of the derivative dis- 
continuity of the XC potential'^^'^'^-s^ 



A. Time Evolution for Quantum Transport 

The time-dependent scheme used in this work is the 
one developed in Ref. 41 and, as in Ref. 63, inter- 
actions in the central region are treated via an ALDA 
from the Bethe Ansatz"^^. For disordered systems, where 
large central regions and configuration averages may be 
needed, such an algorithm may be computationally ex- 
pensive. As described in Sect.VB, a convenient way to 
enhance its numerical efficiency is to use the Lanczos re- 
cursion for time evolution^'^^ (for a quick introduction to 
the Lanczos technique, see Appendix A). 

We start with a concise description of the origi- 
nal algorithm**^, as background to our Lanczos-adapted 
scheme, and we specialize to ID geometries. The nota- 
tion in this and the next section is closer to the one in 
Ref. 41, and thus differs from that in the rest of our 
paper. 

The Hamiltonian we consider is H*°*(i) = H,.; + W(f), 
where W(f) is the external perturbation. In a TDDFT 
approach, the initial, ground state is a single Slater de- 
terminant \^g). It is useful to divide the (ID) space into 
three regions. With s a the site label, we have the region 
L (corresponding to the left lead, with s < — (M -I- 2) ), 
the central region C ( with |s| < M + 1 , i.e the device 
region contains 2M + 3 sites) , and the region R (corre- 
sponding to right lead, with s > [M + 2)). The general 
structure of any bound, extended or resonant one par- 
ticle eigenstate -0 in the Slater determinant \^g) can be 
written as 



-ikis 



+ L_e 



^(s) 



^(s) 



— ikj. 



s < -M-2 

\s\ <M + 1 
s > M + 2 



(25) 



To describe quantum transport, one needs to evolve 
in time the ground state configuration \^g), i.e. each 
one of the single particle eigcnstates ip above. Intro- 
ducing the projection operators Fl,c.r (for example. 



Pl = X^sGL can write (/3 = L,C,R), for the 

generic single particle state, 



|V/3)=P/: 



(26) 



In the same way, we can project the Hamiltonian in the 
difi'erent regions 



H = > Hflfl/, Hflfl'=PflHP«'. 



(27) 



Separating the contribution from the leads in W, the set 
of one-particle equations becomes 

i^\m) = [H(t)+Wws(i)] \m), (28) 

with H(t) = Hoi+Wcc(Oi where Hd is the electron one 
particle Hamiltonian and Wccit) is the external poten- 
tial projected in the central region C. Assuming metallic 
electrodes, 

f 5s,s'WL{t) s<-M-2 
Ww«(i) = < s< |M + 1| . (29) 

[Ss,s'WR{t) s>M + 2 

In the numerical time propagation, the time is dis- 
cretized: tm = 2m5, where 6 is the timestep, m is an 

integer, and the explicit prcfactor 2 is introduced for 
convenience in the formulas. In**^, the one-particle eigcn- 
states are propagated from tm to tm+i using a generalized 
Crank-Nicholson scheme. For the time evolution of each 
one of the one-particle states in I'ig), one gets^^ 



(l + iJH(™)) 



^ ' 2 ^'^ leads 

1 _ iiw^") 



l^(m-H)) 



(l-i(5H(")) 



2 ' ' leads 
^ '2 leads 



1 + 



(m) 
eads 



(30) 



where = \'4'{tm)} and 

H^™) = Hei-h^[Wcc(Wi) + Wcc(im)] (31) 

WL"L = ^[Ww«(Wi) + W;e„da(i™)]. (32) 

Using Eqs.(26,27), and after some algebra, the closed 
equation for the time-evolution in the central region is 



I /("-+!) 



Ic - tSU[J („) 

ic + tsuy 



cx=L,R 



where 



l-zfT4^i"> 
l+i|M/i")' 



3=0 



(33) 

(34) 
(35) 
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and 



.cc-^^ E (36) 



— ±1/- 



The B^"-* matrices have only one non-zero element 

5s.-M--\Ss'-M--i a = L 



- s,s' 



^(0) J Os-M-l^s'-M-l 
Ss,M+lSs,M+l 



a = R 



, (37) 



with 6^°^ = i+Vi+'iS H^ g^j^jj y ^YiQ hopping parameter 
in the leads. The expression for the source state \'y^^) 



and the memory state |Ca"^) are' 



„41. 



Kn)_ 



1 



1 



|Ua), 



\Ua) 



where \ua) is a unit vector such that 

' 1 ^s,M+i a = R 



(38) 
(39) 

(40) 



The scalar quantities Z^"^ and a = L,R are given 
by 



5^ 

j=0 



Wn) _ Z'^, p«a(M+2) , -»z„(M+2)\ y (1 - COS (z^))" 

^ ^ (1 — 2id cos (Zq,)) 



(41) 
(42) 



(1 — 2iS cos (zq)) 



^ (1 - 2i(5cos(za)) 



\ n— J 



and Za = ki for a = L while Za = kr for a = R. For 
n > 2, the quantities 6^"^ in the Eqs. (41,42) are obtained 
by recursion: 



6(0) " 1 + 2(526(0) 
^ :iz,i (6(J) + 6(j-i) + 6(^-2)) 6("-2-i) 

- ^ 22- — 



(43) 



1 + 2525(0) 



and 6("<o) = 0,6(1) ^ i-2gfe|°^ 6(°) and 6(0) the same as 
in Eq.(37). 



B. Lanczos-adapted algorithm 

The basic idea behind the algorithm discussed in the 
previous Section is to discretize the time axis via the 
Crank-Nicholson algorithm before performing the parti- 
tioning in L, C, R regions* 1. One could think of do- 
ing the same using the Lanczos algorithm for the; time 
propagation^o^ (the method is quickly reviewed in Ap- 
pendix A); however, non-commuting parts of the Hamil- 
tonian would appear in the exponent this time, rendering 
formal manipulations more involved. Here, we consider a 
simple shortcut that, while improving the numerical effi- 
ciency of the algorithm of Sect. V A, has the same degree 
of accuracy ( i.e. it is second-order in 6) but avoids work- 



ing with the Lanczos scheme before the partitioningi02_ 
Looking at Eq.(33), we notice that the explicit action of 
Hgg-' occurs in two specific terms: 



IXi) 

1X2) 



ic + isn 



(n) 
eff 



\Ua} 



(44) 
(45) 



'■eff 

where |xi) is the contribution to IV'c*^^^) from the central 
region, and \x2) enters the expressions for the source and 
memory states. For since (5 — >■ 0, one can write, up 
to order two in S 



\xi) = z ^izfylV'c ) 



-ICO- (46) 



For the case of \x2), we define the following quantities: 

1± V3 



which permit to rewrite |x2) as 



1X2) = 



-1 + e 



-iH 



+ e" 



iH 



(n) , 



(47) 



|n„)+0(53) (48) 



If necessary, one can go to higher orders, by imposing 

that (1 + Sx)-'^ = A + T,k ^""'^'^ and finding the coeffi- 
cients A, {a/c} by comparison of the two expressions order 
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by order in 6 (in general, the {ak} will be complex). We 
note that the same Lanczos sequence of basis vectors is 
required for both exponentials in Eq.(48). 
All terms which appear in the propagation scheme of 
Sect. VA and that involve H[,g\ have been re-expressed 
in terms of exponentials, so that Lanczos propagation 
can be used; finally, since H^g'' is complex, Eq.(36), it is 
convenient to split the exponentials; for small S, 

g-2»5H(S) _ e-'5'E„BL'"e-2«*"cce-'''5:„-B™ (49) 

(50) 

For the ID case, the advantage is immediate: the B^"-* 
in Eq.(37) have only one non- vanishing entry and the 
outer exponentials in Eqs. (49,50) reduce to scalars. We 
expect that the splitting will still provide a simplifica- 
tion in the 3D case. To summarize, on increasing the 
size of the central region, our Lanczos adapted scheme 
becomes highly convenient, potentially a significant ad- 
vantage when dealing with disordered and large samples. 



VI. RESULTS 

Recently, a non-perturbative study of finite Anderson- 
Hubbard chains has been performed in terms of Density 
Matrix Renormalization Group and a real-space version 
of Dynamical Mean Field Theory Using different in- 
dicators of delocalization, such as the geometrically av- 
eraged LDOS and the IPR, the main outcome of such 
non-perturbative calculations was a clear indication of a 
tendency to delocalization in a range of U, W values in 
the ground state of these chains. Can a similar behavior 
be observed in shorter chains contacted to semi-infinite 
leads? To address this issue, we present here ground- 
state and dynamical results for short, isolated/contacted 
disordered and interacting chains. No spin-effects are 
considered, i.e. the systems are spin compensated, and 
all applied potentials are spin-independent. The quanti- 
ties we will analyze are the on-site particle density, the 
bond current, and the IPR. In general, we consider box- 
disorder, with onsite energies e; uniformly distributed in 
the interval [—W/2, W/2]. In this case, configuration av- 
erages will be done over a finite number of configurations. 
However, we used binary disorder to discuss the IPR (Sec. 
VIA) and the CPA (Sec. VI B), since, for short chains, 
exact averaging can be done with a manageable num- 
ber of configurations. Further details of each set/type of 
calculations will be provided in the respective sections. 



A. Changing the definition of the IPR 

To analyze our definition of IPR, Eq.(15), we find it 
convenient to consider binary (rather than uniform) dis- 
order in a non-interacting chain of L = 14 sites. For 



IPR 



system 



IPR^ 



Nj_j= 200 



N, =200 




' \ I '^'^ i I '■'"'=-' -0.6 ^ 



0.12 

0.09 

0.06 



-4 -2 



FIG. 1: (Color online) Disorder averaged, energy dependent 
IPR (see main text for the definition) for a Aso-Bso, non- 
interacting binary-alloy chain C with L — 14 sites, contacted 
to in/finite leads. Complete disorder averages were per- 
formed over (^) ~ 3432 configurations. Top-to-bottom pan- 
els show cases with increasing site-disorder W . Left (right) 
panels show results for an IPR estimated on the entire sys- 
tem (only on the chain C), as in Eq.(14) [Eq.(15)]. The 
color/symbol coding in panel a) [b)], also applies to panels 
c),e) [ d),f)]. In the left panels, IPRsystem{ijj) for infinite 
leads {Nid — >■ 00) is not shown. All curves were obtained with 
a small Lorentzian broadening 7, to have a minimum accu- 
racy of 10"* for IPRc{iv) when Nu 00 [black curves m 
panels b),d),f)]. 



binary disorder, choosing L small and even permits i) to 
consider exactly the ^50^50 alloy concentration and ii) 
to perform disorder averages exactly. The chain can be 
isolated or contacted to two ID leads (one at each end of 
the chain); the leads can be finite or semi- infinite (their 
length will be denoted by Ni^). The total number of sites 
in the system is thus M = L + 2Nid. 

We wish at this point to make a short technical di- 
gression on the numerical calculation of the IPR. For the 
energy-dependent IPR, we need the local density of states 
(LDOS) at site(s) i: 



E 



(51) 



where A labels the one-particle eigenstates |A) and eigen- 
values E), of the system. For infinite systems, Eq.(51) 
is not usable directly, and one resorts to Green's func- 
tions [see Eq.(16) in Sect. IVB]. For finite systems, one 
can e.g. use recursion techniques'^, or, as done in this 
Section, perform a direct diagonalization of the disor- 
dered Hamiltonian. However, if the IPR must be deter- 
mined in a range of energies (i.e. for several lo values), 
already for moderate system sizes M , the A- and i-sums 
in Eqs. (14, 15,51) become computationally expensive. In 
Appendix B, we present a technique which permits to 
perform such nested summations in a rather efficient way. 

Results for the IPR according to the two definitions 
C2, Eq.(14), and (3, Eq.(15), are shown in Fig.l (here- 
after, (2 and Cs will be renamed IPRsystem and IPRc, 
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respectively). Calculations with IPRsystem arc reported 
in panels a,c,c). In all three panels, we see that on in- 
creasing the size Nij_ of the (finite) leads, IPRsysiemi'^) 
gets quickly reduced in the region |w| < 2, i.e. for the en- 
ergy range for the extended states in the system (results 
for larger Nid, not shown, confirm this trend). Outside 
the band region, the decay of IPRgystemi^) on increas- 
ing Nid is much slower, and our numerical evidence, to- 
gether with arguments based on the large W limit, shows 
that for larger Nid IPRsystemi^^) vanishes everywhere for 
|a;| > 2 except at the energies of the localized states, 
where it takes the corresponding IPR value. Thus, ir- 
respective of the disorder strength in the finite chain, 
in the limit of semi-infinite leads, IPRsystem indicates 
delocalization for |aj| < 2. However, for a large disor- 
der, the chain becomes disconnected from the leads, and 
this is missed in the vanishing IPRsystem, which simply 
reflect almost everywhere the delocalized states in the 
disconnected leads (as in most quantum transport treat- 
ments, the leads are assumed to be homogeneous and 
non-interacting) . 

In the panels b,d,f) of Fig.l, we show results for 
IPRc{oj): in this case, on increasing Nid, the IPR tends 
to a finite value (the asymptotic value for when Nid ^ 
oo), which better reflects the fact that the localization in 
some region of the system affects the system as a whole 
(we have also verifled that on increasing W, IPRc in- 
creases). Being a local quantity, IPRc obviously de- 
pends on the size and the details of the chain. 

When interactions among particles are taken into ac- 
count, there is another point that is necessary to exam- 
ine. This aspect is specific to our approach to quantum 
transport, where interactions are described within lat- 
tice (TD)DFT. Since the Kohn-Sham system is just a 
fictitious system apt to reproduce the true interacting 
density, a frequency dependent IPR of the Kohn-Sham 
system has actually little physical meaning. Thus, our 
final proposed definition of IPR is: 

iec iec 

Eq.(52), which makes use of the actual particle density, 

can also be used in the interacting and time-dependent 
cases, and thus is both conceptually and operationally 
well defined within (TD)DFT. We have, in the same way 
as for the energy dependent IPR, verified that also this 
definition is sensitive to localization and tends to a finite 
value when the size of the leads tend to infinity. To sum- 
marize, while not arguing that our definition of IPR is 
an optimal or unique indicator of localization in quan- 
tum transport geometries^'', in our simulations we used 
Eq.(52) as a viable prescription. 

B. CPA-DFT for short chains attached to leads 

For a short disordered chain connected to homogeneous 
leads, a treatment based on the single-site CPA amounts 



to introducing a complex, energy- and site-dependent 
self-energy. Numerically, this CPA procedure is quite 
convenient, since it amounts to solving sclf-consistently 
a set of equations. But how accurate is the CPA for 
the kind of (quantum-transport) geometries considered 
in this paper? To answer this question, we have consid- 
ered a disordered chain with L = 14 sites, in absence 
of particle-particle interactions {U = 0). The chain rep- 
resent a A25i?75 and A50-B50 system, i.e., in the nota- 
tion of Section IV B, ca = 0.5 and a ca = 0.25, with 
cb = 1 — Ca- We considered two disorder strengths, i.e. 
W = eA = —£b equal to 0.5 and 1.0. The quantity that 
we intend to examine is the KS average local density of 
states (LDOS), defined as {d{u)))ii = —'!T~^?s{GKs{i^))ii 
{i labels the i-th site in the chain), with {Gks{^) ob- 
tained as described in Section IV B. Strictly speaking, 
the LDOS is not accessible to groimd state DFT; never- 
theless, we think is instructive to look at this quantity 
in the framework of DFT, to compare CPA versus exact 
disorder averaging when interactions are present. For 
numerical convenience, the calculations were performed 
with an additional small Lorentzian broadening. The re- 
sults for [/ = are presented in Fig. 2 (see the figm-e cap- 
tion for additional details). In each panel, the LDOS:s 
are obtained with both exact averages, summing over 
(^.^) — 3432 configurations, and with the CPA. Further- 
more, in each panel, we show LDOS:s at sites adjacent 
to the leads (labeled as " edge" ) and at a site in the cen- 
ter of the chain (labeled as "middle"). Irrespective of 
the strength of the disorder, the results show that, over- 
all, the CPA (at least in this simple single-site formula- 
tion) provides a fair account of the role of the disorder, 
but much of the sharp structures in the exact curves are 
washed out. For example, for W = 1, for w > 1, we 
note a significant depression of the exact LDOS, which 
is completely missed by the CPA. More in general, the 
sharp structures (bound states) outside the continuum 
are removed by the CPA. This can have consequences 
in the long-time limit of qiiantum transport; for pure 
systems, bound-states in the final state KS Hamiltonian 
can give rise to steady-state oscillations^°'^'^°^, but, due 
to disorder as treated in the CPA, such long-lived oscil- 
lations are expected to be damped^°^. The situation is 
less clear for the exactly averaged LDOS: one can ex- 
pect a self- averaging of the current and/or density oscil- 
lations, when the chain increases in size. However, for 
short chains, and for the simple type of binary, on-site 
disorder considered here, (small) long-time oscillations 
could persist. It can be of interest to see at which extent 
this behavior is modified in the presence of inter-particle 
interactions. In Figure 3, the LDOS:s for i7 = and 
U = 3 are plotted for a chain with and without disor- 
der. The curves for [/ = 3 were obtained with lattice 
DFT in the LDA. Starting with the W = case, we 
note that an important the effect of the interactions are 
resonant structures at the top of the band (a; > 2, blue 
curve) . The addition of disorder within a CPA treatment 
has an overall effect similar to what observed in Fig. 2, 
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-1.2 -2.4 -1.6 -0.8 



FIG. 2: (Color online) CPA versus exact averaging for a chain with 14 sites connected to two semi-infinite homogenous leads, 
with transparent boundary conditions. The meaning of each curve is given by the color coding in the top-left panel. In the 
legend, Ex.Av. means complete numerical averaging. The density of states for the leads is nonzero in the interval [—2, 2] and 
the chemical potential jj, — in the leads (half- filling). In each panel, the average LDOS is presented for middle and edge sites 
(see main text), and a obtained with CPA (red and green curves) and exactly (black and orange curves), and two strengths of 
disorder W are considered. For the A50S50 systems (bottom panels) all curves are symmetric w.r.t. u — 0. For W = 1, the 
maxima of some of the sharp peaks at the band-edges are outside the scale shown in the figure. 



namely both the KS sharp structures are dramatically 
smoothened in the KS+CPA LDOS:s (orange curve)^^'^. 
As for the non-interacting case [for reference, the CPA 



w=o, u=o 

W=l, U=0 (CPA) 
W=0, U=3 (KS) 
W=1,U=3(CPA+KS) 
W=l,U=3(Ex. Av. +KS) 




FIG. 3: (Color online) LDOS at an edge site for a disordered 
{A50B50) and interacting 14-site chain contacted to two ho- 
mogenous ID leads at half-filling. The interaction U = 3 and 
the disorder strength W = 1. The other system parameters 
are the same as in Fig. 2. The meaning of each curve is 
given by the color coding in the figure. Also, Ex. Av. means 
complete numerical averaging, while KS indicates that inter- 
actions are treated within a Kohn-Sham scheme. 



result for J7 = is also shown (red curve)], the use of the 
CPA could considerably affect the long-time behavior of 
densities and currents induced by a bias. This is because 
the sharp structures due to the interaction in the DFT- 
LDA (that could induce long time limit oscillations in 
an ALDA treatment) are smoothed by the CPA. For the 
exactly averaged interacting LDOS, we equally observe a 
clear broadening/reduction of the split-off KS structures, 
albeit less pronounced than in the CPA-KS curve. Due 
to the artificial broadening we introduced in our calcu- 
lations, it cannot be excluded that for short chains the 
density and current oscillations would stay long lived. 
The results shown here were obtained from the initial 
state Hamiltonian. However, similar split-off structure 
are present in the case of the final state Hamiltonian, 
and the actual oscillations, independently from the pres- 
ence of disorder, would likely be absent if memory effects 
were taken into consideration^*'". 



C. (TD) DFT results for short chains attached to 
leads: static, transient and steady-state regimes 

In general, for disordered chains attached to semi- 
infinite leads, an exact numerical treatment analogous 
to the one which we will discuss in Sect. VI D is not avail- 
able. In this case, two suitable methods are (TD)DFT 
and the Green's functions technique. These two methods 
are both exact in principle but, in practice, the many- 
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body self-energy in a propagator approach and v^c in 
TDDFT are not known exactly, and approximations arc 
in order. Here we consider a TDDFT description, as pre- 
sented earlier in the paper. Our only (but important) 
approximation will be the use of a local density approx- 
imation (LDA) in the ground state, and its adiabatic 
counterpart (ALDA) during the dynamics. The quality 
of these approximations in the present contexts will be 
briefly discussed at the end of this section. 



1. In equilibrium: The inverse participation ratio. 

We begin studying the system in equilibrium. The 
central region is a chain of L = 15 sites. We consider 
three strengths of disorder in the chain, W = 0,1,3 and 
three values of the interaction, U — 0,2,4. Both U and 
W are given in units of the hopping parameter. In Table 
I, we show the results for the configuration-averaged IPR 
for different densities n = n-^ + in the leads (quarter 
and half- filling cases). The average was performed over 
20 configurations. As discussed in Sect. IV A, the IPR 
is calculated in terms of the KS densities in the central 
region [Eq.(52)]. 

For U = W = Q, the system is homogenous, and we 
have only fully delocalized states and, since L = 15, we 
get an IP Be = L^^ = 0.06. Accordingly, values larger 
than 0.06 would denote a tendency to localization. This 
is what we observe on moving to larger W values, while 
keeping U = 0: localization is maximal for W = 3, both 
at quarter and half filling (however, the degree of localiza- 
tion is different for the two fillings) . A similar dependence 
of the IPR on W is observed for [/ = 2 and !7 = 4. How- 
ever, a different behavior is noted when changing J7 at a 
fixed W (i.e., moving horizontally in the table). Wc see 
that the IPR stays approximatively constant at low W 
but, for larger disorder, the IPR decreases on increasing 
U (see especially the case oiW = 3). That is, increasing 
U decreases localization, a manifestation of the compet- 
ing behavior of interactions and disorder. This has been 
noted before for finite samples (e.g. in terms of exact 
diagonalization^°^ or DMRG calculations^"^). When W 
and U become both very large, calculations as in^^'^'i'^^ 
suggest that localization prevails. From our results, this 
should happen a U values larger than those in the table 
(e.g. for W — 3, a,t U > 4). However, for such inter- 
action strengths, the shortcomings of the Bethe-Ansatz 

TABLE I: Inverse participation ratio in the ground state 





n 


U = 


U = 2 


U = 4 




1.0 


0.06 


0.067 


0.067 


w=o 


0.5 


0.06 


0.067 


0.067 




1.0 


0.068 


0.067 


0.067 


W=l 


0.5 


0.072 


0.071 


0.069 




1.0 


0.077 


0.073 


0.072 


W=3 


0.5 


0.113 


0.096 


0.09 



LDA can become particularly severe^"®. Nevertheless, it 
is quite interesting that a competing regime between dis- 
order and interactions is accounted for within our lattice 
DFT-LDA approach, and with disorder occurring only in 
a sub-region of the system. 

2. Time- dependent densities. 

This behavior should also manifest in the dynamical 
properties of the chain. To see this, we studied the time 
evolution of the system, after the application of a bias in 
the leads. Our choice was to apply the bias only in the 
left lead [i.e., at aU times, 6£)(t) = 0, see Eq.(l)], with 
the following time dependence: 

5^(,) = |Ml + cos.(l+J)]/2 , OJr<T ^53^ 

where T = 3 (all time quantities are expressed in units of 
the reciprocal chain hopping). This choice of T is some- 
what arbitrary, but in this way the effect of 65 (r) is rather 
gradual, a situation expected to be favorable to the use 
of an ALDA based on the Bethe-Ansatz for the ID Hub- 
bard model. Our time-dependent results were obtained 
with propagation time steps of either At = 0.0025 or 
At = 0.0050 and, as in the static case, averaged over 20 
configurations. 

In Fig. 4, we show the effect of disorder and inter- 
actions on the time-dependent density, when bs(t) has 
been applied, with 60 = 1.5. The chemical potential in 
the system was chosen to have half-filling in the leads (i.e. 
n = n-f+n^ = 1). For convenience, we show the deviation 
of the density Arij = ^^(t) — nj(0), rather than the den- 
sity itself, since it illustrates more directly the changes in 
the system. The left panel of Fig. 4 corresponds to when 
neither disorder nor interactions are present in the chain. 
In the transient phase, for sites close to the biased lead, 
we observe a quite sharp rise of the density, while the 
change in density occurs more smoothly for sites closer 
to the unbiasc;d Iciad. It is also clear that the densities in 
the chain attain a steady-state value rather soon, already 
at T < 15, and that there is a quite regular propagation 
of the density front across the chain. Disorder modifies 
in a quite substantially way the situation just described. 
In Fig. 4b), we note an increased An for sites close to 
the biased lead, but the profile of the density propagation 
front is now more irregular and significantly attenuated 
inside the chain. This is also observed for weaker disor- 
der {W = 1, not shown), although the differences from 
the homogeneous case are smaller. 

An interesting fact occurs when introducing interac- 
tions {U = 4) in the chain (Fig. 4c). Now the time- 
dependent density landscape recovers much of the regu- 
larity of the U — W = case, and the propagation of the 
density wave is considerably less attenuated (with respect 
to Fig. 4a, the values of An in the long-time limit are 
reduced). So, it appears that even in the time-dependent 
case, interactions can reduce the effect of disorder. We 
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■a] W=0, U=0 h) W=3, U=0 c) W=3, U=4 




0.0 0.1 0.2 



FIG. 4: (Color online) Difference between non-equilibrium and ground-state (initial) densities, in a chain with L = 15, n = 1 
and box disorder W, when a bias fes(t) is applied in one of the two leads, and bo = 1.5 [see Eq.(53)]. The maximum evolution 
time is Tmax ~ 30, in units of the inverse hopping parameter in the chain. 



already noted such competition of effects when discussing 
the IPR in the ground state, but the results of Fig. 4, 
and other cases we have analyzed, not shown here, con- 
firm the robustness of this behavior with respect to i) 
bias strength (we also considered bg — 0.5), ii) particle 
density (we also investigated the quarter filling regime) 
and, of course, iii) disorder and/or interaction strength. 



u=o 



U=2 



U=4 



3. Time- dependent inverse participation ratio. 



Also accessible within a Kohn-Sham lattice-TDDFT 
scheme is the time-dependent inverse participation ratio, 
which we define via a simple modification of Eq. (52): 



,KS]2 



(r)/E 



(t)]2. Results 



for the time-dependent IPR are shown in Fig. 5 (see the 
figure caption for a definition of all the parameters). 
From Fig. 5, we observe that a larger disorder induces 
a larger IPR, also in the dynamical regime. This holds 
for all cases examined in the figure; at the same time, for 
a fixed disorder strength, interactions make the system 
more delocalized in time (as before, a complete delocal- 
ization corresponds to an IPR ^ 1/L = 0.06. At long 
times, the IPR is reduced compared to its initial value; 
such decrease is almost monotonic for large W, while 
at smaller disorder strengths the IPR grows at first and 
then eventually becomes smaller. The region of increased 
IPR correspond to the transient phase, where the vari- 
ance among the different densities in the chain is largest 
(we have verified that the position at which the IPR at- 
tains its maximum value depends on the way the bias 
is ramped-up). Conversely, the small IPR at long times 
shows that, on average, the densities have the least mu- 
tual variance in a regime where a steady-state current 
can be attained. 
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0.067 
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FIG. 5: (Color online) Time-dependent IPR for a disordered 
and interacting chain with L = 15 sites, contacted to semi- 
infinite leads. Results are shown for different values of box 
disorder (W), interaction ([/), bias (fee) strengths, and for 
different lead densities n. Arithmetical disorder averaging is 
performed over 20 configurations. Panels with one (orange) 
curve refer to the case of W = 3, while those with two (black 
and red) curves refer to W — and W — 1, respectively. 
Panels in the same column pertain to same interaction value 
U, whilst panels in a row refer to common values of bo and 
n, as specified on the right of that row. All panels share the 
same horizontal time interval, but scales on the vertical axes 
are different. 



4- Time- dependent currents and the steady state regime. 

A more direct way to look at the competition of disor- 
der and interactions in the dynamical regime is analyzing 
the behavior of the current (in our case, the charge cur- 
rent). In Fig. 6, we present results for the average current 
at the leftmost-bond in the chain. Altogether, the dif- 
ferent panels of Fig. 6 show the current for several values 
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of disorder/interaction strengths, of the biases, and the 
band fiUing (see the figure caption for details on how the 
results are presented). We start with panel a), corre- 
sponding to [/ — 0, a bias bo — 0.5, and a density (in 
the leads) n = 0.5 (quarter filhng). For W — (black 
curve), the current reaches its steady-state value after a 
relatively short transient. On adding disorder {W = 1, 
red curve), the length of transient increases, but eventu- 
ally a steady state is reached. On further increasing W, a 
steady-state current is not reached within the simulation 
window (orange curve). However, the current in Fig. 6 
(orange curve) and the currents at the other bonds in 
the chain (not-shown) become progressively close to the 
same average value (with smooth and decaying oscilla- 
tions), suggesting that a steady state is reached outside 
our simulation window. The overall trend in the panel 
a) is that an increasing disorder reduces the long-time 
(steady-state) value of the current. Analogous behavior 
is seen in panels a), b) of the U — column, which corre- 
spond to different choices of the bias bo and the density n 
in the leads. Moving to the other two columns {U = 2,4), 
we see that the changes in each panel follow the same 
pattern, namely when disorder increases the current gets 
reduced. 

A different perspective emerges from Fig. 6 when we 
mutually compare panels within the same row [e.g. pan- 
els a), d), c)]. In this case, for a given value of disorder 
strength, bias and density, the current increases at larger 
values of U, an effect of competing disorder and interac- 
tion in the non-equilibrium regime. In the first row of 



U=0 U=2 U=4 




time 

FIG. 6: (Color online) Time-dependent average currents for 
a chain with L = 15 sites, and different interactions strengths 
{U = 0,2,4). The currents shown are computed at the left- 
most bond in the chain. On the horizontal axis, the time is 
expressed in units of the inverse hopping parameter. The bias 
and the band filling for panels in the same row are specified 
on the right, while in each panel, the current is displayed for 
three disorder strengths {W = 0,1,3). Color coding for all 
panels is specified in the right lowermost panel. The con- 
figuration averages were obtained from 20 instances of box 
disorder. 



Fig. 6 the current is always increasing when U becomes 
larger, independently of the value of W. Conversely, in 
the other two rows, depending on the value of W, the cur- 
rent can have a non-monotonic dependence on U. Such 
non-monotonic (dynamical) behavior is consistent with 
results from ground-state studies (see e.g.^"''), and its 
dependence on quantities such as the density n or the 
disorder strength W is plausible. To have a more com- 
plete picture of this tendency, calculations for larger U 
values should be performed; however, as mentioned ear- 
lier, for those U values, the utility of an ALDA-TDDFT 
approach would be significantly diminished. 

So far, we have not mentioned at all the oscillatory 
behavior of the current in some of the panels of Fig. 6. 
Current and/or density oscillations in the long-time limit 
can result from different factors, such as single-particle 
bound states^^^'^*'^ (in our case, due to diagonal disor- 
der), discontinuities in the XC potential^'^, or a slosh- 
ing motion of the charge between different regions of the 
device and/or the terminal sites of the leads. This lat- 
ter mechanism has been pointed out and analyzed in a 
study^^^ of the dynamical effects of image charge in quan- 
tum transport. In our system, we have observed that dif- 
ferent and independent sets of oscillations can emerge in 
different parts of the chain, due to the inhomogeneities 
in the energetic landscape introduced by disorder (cur- 
rent and cumulative density oscillations are more/less 
pronounced in different subregions of the chain). In re- 
cent work^odoo.iii^ Hartree-Fock and TDDFT-ALDA ap- 
proaches have been compared to results from Kadanoff- 
Baym dynamics. These studies clearly point out the im- 
portance of nonlocal (in space and time) contribution be- 
yond the instantaneous density, and suggest^''^'^^^ that 
current /density oscillations of the kind mentioned above 
are likely fragile against the inclusion of memory effects. 
We wish to add here that, for disordered systems, an- 
other obvious cause of dephasing of the oscillations is 
disorder averaging. This can be already argued at the 
non-interacting level. While bound-states can certainly 
be present in a specific instance of disorder, the induced 
oscillations are most probably to be washed out by config- 
uration averaging, due to destructive interference among 
the oscillation from different configurations. This re- 
ceives indirect support from observing that already in the 
ground state of a device with binary disorder (Fig. 3), 
split-off structures are largely reduced or washed out by 
exact or CPA averaging. In our time-dependent calcula- 
tions, as commonly done in the literature, configuration 
averages are based on a limited (and thus incomplete) 
number of random disorder realization; in this case the 
aforementioned cancellation effects of the oscillations can 
be incomplete. In any case, the feature emerging from 
our calculations, namely a competition between disorder 
and interaction, appears to unrelated to current oscilla- 
tions, since it present itself also when such oscillations 
are missing [see e.g., in panels a), d), g) in Fig. 6]. 

As a final, but important remark, we observe that our 
results for the current lack of "reciprocity", i.e. for a 
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fixed U, and the W:s considered, the current is mono- 
tonically decreasing as a function of W , while a com- 
peting behavior could be expected from studies in the 
ground state^"'*. This can possibly be due to a limitation 
of the ALDA (and not of TDDFT, which is in principle 
an exact theory). At the same time, we wish to point out 
that, out of equilibrium and in the non-linear regime, the 
strength of the bias 60 can also significantly influence the 
competition between disorder and interactions (results 
for 60 = 1.5, n = 0.5, not shown here, compared to those 
in panels a),d),g), are consistent with this assertion). 

D. Isolated short chains: TDDFT and exact 
results. 




time time 



FIG. 7: (Color online) Exact (left panels) vs lattice TDDFT 
(right panels) averaged density results for an isolated chain 
with i = 18 sites and U = 2. The chain has a central region of 
8 sites, with disorder and interaction, and two finite "leads" of 
5 sites each (see the schematic rendering). In the chain there 
are N = A, spin-compensated particles and the dynamics is 
induced with a time dependent "bias" which is uniformly ap- 
plied to the nine-leftmost sites. The time dependence of the 
bias is the same of Eq.(53), with strength 6o=0.5. The color 
coding in the top-left panel applies to all panels in the fig- 
ure. The disorder averaged densities were obtained using 100 
disorder configurations. 

The quantum transport results just discussed are based 
on approximate XC potentials, obtained via an (adia- 
batic) local density approximation to v^c- Comparisons 
between exact and (A)LDA results have been performed 
before in different contexts'^^'^^'^^'^°^^*'°, but without con- 
sidering the case of interest here, namely configuration 
averaged quantities in disordered and interacting sam- 
ples contacted to semi-infinite leads. 

An ideal way of doing this would be using time- 
dependent DMRG results as benchmark, to obtain a com- 
prehensive assessment of lattice (TD)DFT for the tran- 
sient behavior of current and densities. From the com- 



putational point of view, this kind of investigations are 
expected to be rather expensive, and we are not aware of 
any published work on the subject. Here we take a much 
simpler view, and consider small isolated disordered and 
isolated chains with few electrons, which can be treated 
via exact diagonalization. An extensive comparison be- 
tween exact and lattice (TD)DFT results is outside the 
scope of this work, and we briefiy discuss just one exam- 
ple, using a single case to gain some general insight. 

In Fig. 7, we examine the time-dependent densities, for 
a short isolated cluster (all the parameters are specified 
in the figure caption). We mention that, differently from 
the previous sections, the "bias" is applied on the left half 
of the cluster, i.e. also in part of the interacting, disor- 
dered region (see the schematic rendering in at the top of 
Fig. 7). Also, due to the small size of our system, several 
reflections at the cluster boundary occur during the simu- 
lation interval considered in Fig. 7. Without going in any 
detail, the salient features emerging from Fig. 7 are that 
TDDFT performs better at larger disorder strength, and 
that, compared to the exact results, TDDFT enhances 
the spread among the time-evolved densities. From this 
example, it is apparent that lattice (TD)DFT is able to 
reproduce the qualitative behavior of the exact results, 
but good quantitative agreement is lacking. More inves- 
tigations in this direction and an improvement of the XC 
potentials are certainly required. 



VII. CONCLUSIONS 

We have used lattice TDDFT to study the quantum 
transport properties of short, disordered and interact- 
ing chains contacted to semi-infinite leads. Our work is 
largely exploratory in character, since we have addressed 
only superficially several issues connected to the non- 
equilibrium physics of disordered interacting systems. 

In principle TDDFT is an exact approach, but (at 
times severe) approximations are usually made for the 
key quantity of the approach, the exchange-correlation 
potential. We have employed one of them, the adia- 
batic local density approximation (ALDA). Disordered 
systems, with a strongly varying local environment, are 
quite difficult tests for the ALDA which, however, by 
comparison to exact benchmarks, certainly appears ap- 
propriate at the qualitative level for not too-fast varying 
time-dependent perturbations. 

Within these boundaries, we have been able to address 
a type of system which is not easily accessible, for a rea- 
son or another (limitations in principle, numerical costs, 
dimensionality, etc.) to several of the other methods cur- 
rently available. In fact, we are not aware of any existing 
work for the quantum transport geometries considered 
here, where, in the presence of semi-infinite contacts, cur- 
rents and densities in the disordered and correlated sam- 
ple have been followed in time from the initial transient 
phase to the long-time (possibly steady-state) regime. To 
perform our study, we have introduced some modifica- 
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tions to the formalism, and modified the definition of 
one of the standard indicators of localization, the inverse 
participation ratio. We have also explored the perfor- 
mance of the coherent potential approximation, one of 
the popular schemes to perform disorder averaging. The 
sharp spectral fciatures due to disorder in the device are 
exceedingly smoothed by the CPA, and this can affect 
the behavior of the averaged steady-state currents. 

Our time-dependent results show rather evident signa- 
tures of dynamical delocalization due to the dynamical 
competition of disorder and interaction in the sample. 
This is consistent with the qualitative picture usually 
adopted for systems in the ground state, where inter- 
actions produce a "screening" of the disorder - i.e., the 
"attractive" behavior of low-energy impurities is com- 
pensated by the local repulsive interactions, thus pro- 
viding a less corrugated energy landscape. However, in 
the present case, an additional role is played by the (non- 
weak) electric bias in the lead(s), which significantly (and 
dynamically) modifies the ground state energy landscape 
and wavef unctions. At which extent our findings remain 
robust against an improvement of the XC potentials is at 
present an open question. This certainly calls for better 
potentials and future investigations^^^, hopefully com- 
bined also with other methods. 

We thank C.-O. Almbladh for valuable discussions 
and for critically reading the manuscript. This work 
was supported by European Theoretical Spectroscopy 
Facility (INFRA-2007-211956). 



APPENDIX A: Lanczos time propagation 

We briefly summarize the Lanczos method, as given 
in Ref. 103. A useful comparative study between the 
Lanczos method and other integration schemes can be 
found in Ref. 113. Consider a system described by a 
TD Hamiltonian H{t). If, for example, we use the mid- 
point approximation for the time propagator and wish 
to evolve the system in the time interval {t + A,t), we 
obtain 



(54) 



where is the (known) initial wavefunction. Consider 
a finite Lanczos sequence {114}}, obtained by starting 
acting on the 'seed' |$() = |Vb). Using {|Vfe)} as a trun- 
cated basis, we get 



Ml 



(55) 



k=0 



where is the tridiagonal representation for H{t+ A/2) 
in such a basis. Inserting a complete set of eigenstates 
for the truncated space, Hl\X) = ex\X), 



K 



fe=0 



(A|Vo) 



(56) 



where |^'t+A) is finally expressed in the basis of the orig- 
inal many body Hamiltonian. The method requires a 
partial orthogonalization on the fly of the Lanczos basis 
in order to preserve accuracy along the trajectory. For 
a simple estimate of the truncation error in Eq.(55), see 
the discussion in Ref. 103. 



APPENDIX B: Computationally efficient Lorentzian 



To determine the energy dependent IPR of a finite sys- 
tem using the Nx eigenstatcs/values of the Hamiltonian, 
we have to determine the local density of states (LDOS) 
n{(jj) = e^5{uj — ex). This has to be done for the 
values of the chosen energy grid. Furthermore, in some 
cases (as in our work here) the IPR is also averaged over 
Nd disorder configurations. Since the system is finite, it 
is expedient to introduce a Lorentzian broadening F, and 
to define a broadened LDOS 



nF{u)) = n{uj) * ^'"(w) 



A 



c^L'^iuj-ex), (57) 



where L^{uj) = (r/7r)(a;2 + F^)"!. When Nx,N^,Nd 
are large, to compute •nF{u)) directly from Eq.(57) can 
be computationally intensive. Fast Fourier Transform 
(FFT) is the method of choice in these cases, but it 
requires uniform sampling, whilst the eigenvalues poles 
ex are in general unevenly distributed. This issue can 
be avoided with the approach described here, which is 
fast and accurate, and in principle should be relevant 
to FFT integration problems for large data set from 
nonuniform, adaptive or curvilinear sampling^ To 
begin with, it is convenient to have a smooth function 
C^{u)) that decays rapidly in both u) and t spaces, min- 
imizes ui sampling and allows smooth uniform sampling 
of 'nP^ (uj) = n{u) * C^iuS). Then, the actual procedure 
is as follows: i) we sample nP"" {uj) ^ {ui) , C~'{u})-, ii) 
via FFT, we compute n'^'' (t) , (t) , C^{t)-, in) we obtain 
n^[t) = [n'='''{t)/C''{t)] X L^{t); iv) via inverse FFT, we 
compute n^(w). 

The softening function should be positive definite in 
both ui, t spaces, to ensure easy deconvolution (we discard 
band limited functions such as rectangle, triangle, etc., 
since they require careful location of their zeroes). We 
also discard C' = cxp(— 70;^), since a) C''{t) may decay 
too fast, with deconvolution instabilities where L^{t) is 
still non negligible; b) it is still expensive to determine 
n(w) *exp(— 7w^) on samphng points. Based on these 
considerations, our optimal choice is 



6-^(0;) = ^(1-7 



d_ 
(?7 



)e-'^l"l (X e"''l"l *e" 



-7|w| 



(58) 



The function C^{ij) decays exponentially, i.e. is "prac- 
tically" band limited with a small sampling domain; 
C'^(t) oc (f^ -f- 7^)~^ is always positive, and it decays 
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slower than L^{t), thus avoiding deconvolution instabil- 
ities. Also, as results of the self-convolution in the def- 
inition of C', the cusp e~'''l'^~^^l, when ex is off-grid, is 
smoothed, reducing sampling errors. Finally, as a crucial 
advantage of the method, n*-^ (w) can be computed re- 
cursively, needing the A-sum just once (instead of at all 
the sampling points). In fact, writing n'^'' (uj) = n+{ijj) + 
with n±(w) = '£xC^(^^'^^'^~'^^^iH'^\ - we 



get (9 is the Heaviside function): 

h±{uj T A) = e~T"^[n±(w) + 7An±(a;)] 
+(1 - 7^) X y c^e±^('^-^^T^5 

X &{±{oj-ex))e{A±{ex-oj)). (59) 
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